package fem2;

import inf.jlinalg.IVector;
import inf.jlinalg.SolveFailedException;
import inf.text.ArrayFormat;
import junit.framework.Assert;
import math2.FunctionRnToR;
import math2.RampFunctionRToR;

import org.junit.Test;

import fem2.analysis.Analysis;
import fem2.analysis.LinearStaticAnalysis;
import fem2.element.StructuralElement;
import fem2.enu.EchoLevelType;
import fem2.enu.MeshPartType;
import fem2.enu.State;
import fem2.material.StVenantKirchhoffMM;

public class PatchTestQ8 {

	public Model createModel() {
		/*
		 * mesh
		 */
		Mesh mesh = new Mesh();

		double l = 1;
		double h = 1;
		Node n1 = mesh.addNode(new Node(0, 0));
		Node n2 = mesh.addNode(new Node(l, 0));
		Node n3 = mesh.addNode(new Node(l, h));
		Node n4 = mesh.addNode(new Node(0, h));
		Node n5 = mesh.addNode(new Node(l / 2, 0));
		Node n6 = mesh.addNode(new Node(l, h / 2));
		Node n8 = mesh.addNode(new Node(l / 2, h));
		Node n9 = mesh.addNode(new Node(0, h / 2));

		Face f1 = mesh.addFace(new Face(n1, n2, n3, n4, n5, n6, n8, n9));

		mesh.addEdge(new Edge(n2, n3, n6));

		mesh.setNumberOfFields(2);
		/*
		 * element and model
		 */
		double E = 2;
		double nu = 0.3;
		double t = 1.0;
		State ss = State.PLANE_STRESS;
		MaterialModel cm = new StVenantKirchhoffMM(E, nu, t, 0, ss);

		StructuralElement e1 = new StructuralElement(f1, cm);
		// e1.setIntegrationType(IntegrationType.REDUCED_INTEGRATION);

		Model m = new Model(mesh);
		m.addElement(e1);

		/*
		 * constraint
		 */
		Constraint c1 = new Constraint(false, false);
		Constraint c2 = new Constraint(true, false);
		Constraint c3 = new Constraint(false, true);
		n1.setConstraint(c1);
		n2.setConstraint(c2);
		n5.setConstraint(c2);
		n4.setConstraint(c3);
		n9.setConstraint(c3);

		return m;
	}

	@Test
	public void testEnumeration() {
		Model m = new PatchTestQ8().createModel();
		m.preAnalysis();
		System.out.println(ArrayFormat.format(m.getMesh().getFace(0).getDofNumbers()));
	}

	// @Test
	// public void testKe() {
	// Model m = new PatchTestQ8().createModel();
	// m.preAnalysis();
	// System.out.println(MatrixFormat.format(m.getElement(0).computeKte()));
	// }

	@Test
	public void testGaussPoint() {
		MeshPartType type = MeshPartType.QUAD4;

		int p = 2;
		int ng = IntegrationManager.getInstance().getNumIntegrationPoints(type, p);

		for (int i = 0; i < ng; i++) {
			double[] xi = IntegrationManager.getInstance().getIntegrationPoint(type, p, i);
			double w = IntegrationManager.getInstance().getIntegrationWeight(type, p, i);
			Debugger.watch(xi);
			Debugger.watch(w);
		}
	}

	double[][] Uh = new double[][] { { 0, 0 }, { 0.5, 0 }, { 0.5, -0.15 }, { 0, -0.15 },
			{ 0.25, 0 }, { 0.5, -0.075 }, { 0.25, -0.15 }, { 0, -0.075 } };

	@Test
	public void testSolve1() throws SolveFailedException {
		Model m = new PatchTestQ8().createModel();
		double F = 0.5;
		Force l1 = new Force(m.getMesh().getNode(1), F / 3, 0);
		Force l2 = new Force(m.getMesh().getNode(5), 4 * F / 3, 0);
		Force l3 = new Force(m.getMesh().getNode(2), F / 3, 0);
		m.addLoad(l1);
		m.addLoad(l2);
		m.addLoad(l3);

		Analysis an = new LinearStaticAnalysis(m);
		an.getStrategy().setEchoLevel(EchoLevelType.OUTPUT);
		an.run();

		/*
		 * print the displacement and check with theoretical values
		 */
		for (int i = 0; i < m.getMesh().countNodes(); i++) {
			double[] uHatLocal = m.getMesh().getNode(i).getUHatLocal();
			System.out.println("node " + (i + 1) + ":" + ArrayFormat.format(uHatLocal));
			Assert.assertEquals(Uh[i][0], uHatLocal[0], 1e-10);
			Assert.assertEquals(Uh[i][1], uHatLocal[1], 1e-10);
		}

		for (int i = 0; i < m.countElements(); i++) {
			StructuralElement e = (StructuralElement) m.getElement(i);
			for (int j = 0; j < e.getNumIntegrationPoints(); j++) {
				IVector sigma = e.getStress(j);
				Assert.assertEquals(1.0, sigma.get(0), 1e-10);
				Assert.assertEquals(0.0, sigma.get(1), 1e-10);
				Assert.assertEquals(0.0, sigma.get(2), 1e-10);
			}
		}
	}

	@Test
	public void testSolve2() throws SolveFailedException {
		Model m = new PatchTestQ8().createModel();
		double F = 1;
		Mesh mesh = m.getMesh();

		FunctionRnToR ld = new RampFunctionRToR(-1, 1, F, F);
		Load el = new DistributedLoad(mesh.getEdge(0), ld, null);
		m.addLoad(el);

		Analysis an = new LinearStaticAnalysis(m);
		an.getStrategy().setEchoLevel(EchoLevelType.OUTPUT);
		an.run();

		/*
		 * print the displacement and check with theoretical values
		 */
		for (int i = 0; i < m.getMesh().countNodes(); i++) {
			double[] uHatLocal = m.getMesh().getNode(i).getUHatLocal();
			System.out.println("node " + (i + 1) + ":" + ArrayFormat.format(uHatLocal));
			Assert.assertEquals(Uh[i][0], uHatLocal[0], 1e-10);
			Assert.assertEquals(Uh[i][1], uHatLocal[1], 1e-10);
		}

		for (int i = 0; i < m.countElements(); i++) {
			StructuralElement e = (StructuralElement) m.getElement(i);
			for (int j = 0; j < e.getNumIntegrationPoints(); j++) {
				IVector sigma = e.getStress(j);
				Assert.assertEquals(1.0, sigma.get(0), 1e-10);
				Assert.assertEquals(0.0, sigma.get(1), 1e-10);
				Assert.assertEquals(0.0, sigma.get(2), 1e-10);
			}
		}
	}

	@Test
	public void testSolve3() throws SolveFailedException {
		Model m = new PatchTestQ8().createModel();
		Mesh mesh = m.getMesh();

		double du = 0.5;
		Constraint c1 = new Constraint(false, false);
		c1.setValue(0, du);
		mesh.getNode(1).setConstraint(c1);
		Constraint c2 = new Constraint(false, true);
		c2.setValue(0, du);
		mesh.getNode(2).setConstraint(c2);
		mesh.getNode(5).setConstraint(c2);

		Analysis an = new LinearStaticAnalysis(m);
		an.getStrategy().setEchoLevel(EchoLevelType.OUTPUT);
		an.run();

		/*
		 * print the displacement and check with theoretical values
		 */
		for (int i = 0; i < m.getMesh().countNodes(); i++) {
			double[] uHatLocal = m.getMesh().getNode(i).getUHatLocal();
			System.out.println("node " + (i + 1) + ":" + ArrayFormat.format(uHatLocal));
			Assert.assertEquals(Uh[i][0], uHatLocal[0], 1e-10);
			Assert.assertEquals(Uh[i][1], uHatLocal[1], 1e-10);
		}

		for (int i = 0; i < m.countElements(); i++) {
			StructuralElement e = (StructuralElement) m.getElement(i);
			for (int j = 0; j < e.getNumIntegrationPoints(); j++) {
				IVector sigma = e.getStress(j);
				Assert.assertEquals(1.0, sigma.get(0), 1e-10);
				Assert.assertEquals(0.0, sigma.get(1), 1e-10);
				Assert.assertEquals(0.0, sigma.get(2), 1e-10);
			}
		}
	}

	public static void main(String[] args) throws SolveFailedException {
	}
}
